bosaFS analyses

## 
##         appendicularia    calanus glacialis c           calanus sp n 
##             0.01610809             0.02767183             0.18542563 
##       cyclopoidae sp c                    egg                  other 
##             0.13284103             0.18323460             0.05019070 
##    other calanoid sp c     other calanus sp c other copepodite and n 
##             0.02466932             0.02191025             0.04637669 
##     pseudocalanus sp c     pseudocalanus sp n 
##             0.08650491             0.22506695

##         appendicularia    calanus glacialis c           calanus sp n 
##             0.05119706             0.32603872             0.07813962 
##       cyclopoidae sp c                    egg                  other 
##             0.02777306             0.02417236             0.02319161 
##    other calanoid sp c     other calanus sp c other copepodite and n 
##             0.06510235             0.30352209             0.02277575 
##     pseudocalanus sp c     pseudocalanus sp n 
##             0.06348159             0.01460578

#Beaufort Sea: Mackenzie Shelf, Amundsen Gulf Maud
#Kitikmeot: Coronation Maud, Larsen Sound, Peel Sound
#Baffin Bay: North Water, Lancaster Sound, NW Baffin Bay
dfish$province <- with(dfish, ifelse(region %in% c("Mackenzie Shelf", "Amundsen Gulf Mouth"), "Beaufort Sea",
                                    ifelse(region %in% c("Coronation Maud", "Larsen Sound - Victoria Strait", "Peel Sound"), "Kitikmeot",
                                           ifelse(region == "NEG", "Greenland Sea",
                                                  "NW Baffin Bay"))))
reg_initials <- data.frame(region = levels(dfish$region), initials = c("AGM", "CM", "LS", "LV", "MS", "NEG", "NW", "PS", "WBB"))

dfish <- merge(dfish, reg_initials)
gutcontent2 <- merge(gutcontent2, reg_initials)

dfish$region_year <- as.factor(paste(dfish$initials, dfish$year))
gutcontent2$region_year <- as.factor(paste(gutcontent2$initials, gutcontent2$year))

prey_matrix_comb2 <- data.table::dcast(data.table::setDT(gutcontent2), unique_fish_id~prey_id_comb, length) %>% merge(., dfish[,c("unique_fish_id", "region_year", "region", "station", "open_water_day","surf_sal_kgm3", "surf_temp_degC", "NASC_zoo", "indstandard", "prof_mel", "est_standard_length", "fish_cond", "feeding_success")]) %>% select(-unique_fish_id)

dfish$region_year <- ordered(dfish$region_year, levels = c("MS 2009","MS 2010","MS 2014","MS 2015",
                                                           "AGM 2015", 
                                                           "CM 2011" ,"CM 2016","CM 2017","CM 2018",
                                                           "LV 2010", "LV 2016","LV 2018",
                                                           "PS 2010","PS 2011","PS 2014","PS 2015","PS 2016",
                                                           "LS 2010","LS 2011","LS 2015" ,"LS 2016","LS 2017",
                                                           "NW 2014","NW 2016", "NW 2018",
                                                           "WBB 2015","WBB 2016",
                                                           "NEG 2017"
                                                           ))
station_data <- aggregate(data = dfish, unique_fish_id ~ region_year + station + surf_temp_degC + surf_sal_kgm3 + open_water_day + NASC_zoo + indstandard + prof_mel, length)

Section 1. Analyses multivariées

Environnement

PCA des régions années et conditions physiques

env_pca <- station_data%>% select(surf_temp_degC, surf_sal_kgm3, open_water_day, region_year, prof_mel) %>% PCA(quali.sup = c(4), graph = FALSE)

station_data$pc1 <- env_pca$ind$coord[, 1]
station_data$pc2 <- env_pca$ind$coord[, 2]

pca.vars <- env_pca$var$coord %>% data.frame
pca.vars$vars <- rownames(pca.vars)

posArrowsX <- c(4, 4, 4,4)
posArrowsY<- c(4,4,4,4)
colours <- c(pals::stepped(n=24)[13:16], pals::polychrome(n=4)[4], pals::stepped(n=24)[5:11], pals::stepped2(n=20)[13], pals::stepped(n=24)[c(1:4,17:20)],pals::stepped3(n=16)[16],pals::stepped(n=24)[21:23], pals::stepped3(n=20)[c(5:6)], "red")
names(colours) <- levels(station_data$region_year)
gg_PCA_env <- ggplot() +  
  geom_point(data = station_data, aes(x = pc1, y = pc2, colour = region_year), size = 2.5)+
  #scale_colour_manual(values = colours)+
  coord_equal()+
  theme_classic()+
  theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
        plot.background = element_rect(fill = "transparent", color = NA),
        legend.background = element_rect(fill=alpha('white', 0.4)))+
  labs(x= paste("PC1 (", round(env_pca$eig["comp 1", "percentage of variance"], digits = 2), "%)"),
       y= paste("PC2 (", round(env_pca$eig["comp 2", "percentage of variance"], digits = 2), "%)")) +
  geom_hline(yintercept = 0, lty = 2) +
  geom_vline(xintercept = 0, lty = 2) +
  geom_segment(data = pca.vars, aes(x = 0, xend = Dim.1*3, y = 0, yend = Dim.2*3),
               arrow = arrow(length = unit(0.025, "npc"), type = "open"), 
               lwd = 0.5) + 
  geom_text(data = pca.vars, 
            aes(x = Dim.1*posArrowsX, y =  Dim.2*posArrowsY, label = vars), 
            check_overlap = FALSE)+
  scale_color_manual(values = colours)
gg_PCA_env  

Les clusters se font par région-année. L’Ouest et l’est semblent se mélanger, mais l’Arctique canadien (Kitikmeot) semble relativement distinct dans ses caractéristiques océanographiques. Le Mackenzie Shelf en 2010 se démarque particulièrement dû aux fortes températoires et sa faible profondeur de mélange. 2 stations de Peel Sound en 2014 semblent être particulièrement froides avec une couche de mélange très profonde, tout comme une statioin du NEG. Les deux stations du Lancaster Sound en 2010 se démarquent du cluster de LS.

PCA des régions-années, conditions environnementales et zoo

envzoo_pca <- station_data%>% select(surf_temp_degC, surf_sal_kgm3, open_water_day, NASC_zoo, region_year, prof_mel) %>% PCA(quali.sup = c(5), graph = FALSE)

station_data$pc1_zoo <- envzoo_pca$ind$coord[, 1]
station_data$pc2_zoo <- envzoo_pca$ind$coord[, 2]

pca.vars <- envzoo_pca$var$coord %>% data.frame
pca.vars$vars <- rownames(pca.vars)

posArrowsX <- c(4, 4, 4,4,4)
posArrowsY<- c(4,4,4,4,4)
gg_PCA_envzoo <- ggplot() +  
  geom_point(data = station_data, aes(x = pc1_zoo, y = pc2_zoo, colour = region_year), size = 2.5)+
  #scale_colour_manual(values = colours)+
  coord_equal()+
  theme_classic()+
  theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
        plot.background = element_rect(fill = "transparent", color = NA),
        legend.background = element_rect(fill=alpha('white', 0.4)))+
  labs(x= paste("PC1 (", round(envzoo_pca$eig["comp 1", "percentage of variance"], digits = 2), "%)"),
       y= paste("PC2 (", round(envzoo_pca$eig["comp 2", "percentage of variance"], digits = 2), "%)")) +
  geom_hline(yintercept = 0, lty = 2) +
  geom_vline(xintercept = 0, lty = 2) +
  geom_segment(data = pca.vars, aes(x = 0, xend = Dim.1*3, y = 0, yend = Dim.2*3),
               arrow = arrow(length = unit(0.025, "npc"), type = "open"), 
               lwd = 0.5) + 
  geom_text(data = pca.vars, 
            aes(x = Dim.1*posArrowsX, y =  Dim.2*posArrowsY, label = vars), 
            check_overlap = FALSE)+
  scale_color_manual(values = colours)
gg_PCA_envzoo

On a encore MS 2010 qui se distingue particulièrement. L’Arctique canadien se discrimine un peu moins, malgré que CM, LV et PS semblent appartenir au même groupe.

PCA des régions-années, conditions environnementales et données de filets de larves

envfish_pca <- station_data%>% select(surf_temp_degC, surf_sal_kgm3, open_water_day, indstandard, region_year, prof_mel) %>% PCA(quali.sup = c(5), graph = FALSE)

station_data$pc1_fish <- envfish_pca$ind$coord[, 1]
station_data$pc2_fish <- envfish_pca$ind$coord[, 2]

pca.vars <- envfish_pca$var$coord %>% data.frame
pca.vars$vars <- rownames(pca.vars)

posArrowsX <- c(4, 4, 4,4,4)
posArrowsY<- c(4,4,4,4,4)
gg_PCA_envfish <- ggplot() +  
  geom_point(data = station_data, aes(x = pc1_fish, y = pc2_fish, colour = region_year), size = 2.5)+
  #scale_colour_manual(values = colours)+
  coord_equal()+
  theme_classic()+
  theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
        plot.background = element_rect(fill = "transparent", color = NA),
        legend.background = element_rect(fill=alpha('white', 0.4)))+
  labs(x= paste("PC1 (", round(envfish_pca$eig["comp 1", "percentage of variance"], digits = 2), "%)"),
       y= paste("PC2 (", round(envfish_pca$eig["comp 2", "percentage of variance"], digits = 2), "%)")) +
  geom_hline(yintercept = 0, lty = 2) +
  geom_vline(xintercept = 0, lty = 2) +
  geom_segment(data = pca.vars, aes(x = 0, xend = Dim.1*3, y = 0, yend = Dim.2*3),
               arrow = arrow(length = unit(0.025, "npc"), type = "open"), 
               lwd = 0.5) + 
  geom_text(data = pca.vars, 
            aes(x = Dim.1*posArrowsX, y =  Dim.2*posArrowsY, label = vars), 
            check_overlap = FALSE)+
  scale_color_manual(values = colours)
gg_PCA_envfish

PCA des régions-années, et relations entre l’environnement et l’abondance de zoo et de larves

envbio_pca <- station_data%>% select(surf_temp_degC, surf_sal_kgm3, open_water_day, NASC_zoo,indstandard, region_year, prof_mel) %>% PCA(quali.sup = c(6), graph = FALSE)

station_data$pc1_bio <- envbio_pca$ind$coord[, 1]
station_data$pc2_bio <- envbio_pca$ind$coord[, 2]

pca.vars <- envbio_pca$var$coord %>% data.frame
pca.vars$vars <- rownames(pca.vars)

posArrowsX <- c(4, 4, 4,4,4,4)
posArrowsY<- c(4,4,4,4,4,4)
gg_PCA_bio <- ggplot() +  
  geom_point(data = station_data, aes(x = pc1_bio, y = pc2_bio, colour = region_year), size = 2.5)+
  #scale_colour_manual(values = colours)+
  coord_equal()+
  theme_classic()+
  theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
        plot.background = element_rect(fill = "transparent", color = NA),
        legend.background = element_rect(fill=alpha('white', 0.4)))+
  labs(x= paste("PC1 (", round(envbio_pca$eig["comp 1", "percentage of variance"], digits = 2), "%)"),
       y= paste("PC2 (", round(envbio_pca$eig["comp 2", "percentage of variance"], digits = 2), "%)")) +
  geom_hline(yintercept = 0, lty = 2) +
  geom_vline(xintercept = 0, lty = 2) +
  geom_segment(data = pca.vars, aes(x = 0, xend = Dim.1*3, y = 0, yend = Dim.2*3),
               arrow = arrow(length = unit(0.025, "npc"), type = "open"), 
               lwd = 0.5) + 
  geom_text(data = pca.vars, 
            aes(x = Dim.1*posArrowsX, y =  Dim.2*posArrowsY, label = vars), 
            check_overlap = FALSE)+
  scale_color_manual(values = colours)
gg_PCA_bio

On a encore la démarcation Kitikmeot reste du groupe, aninsi que MS 2010 qui est seul. L’est et l’ouest se ressemblent beaucoup avec quelques stations qui font exception.

Ajoutons les larves

fish_pca <-  dfish %>% select(region_year, feeding_success, surf_temp_degC, surf_sal_kgm3, open_water_day, NASC_zoo,indstandard, prof_mel, est_standard_length, fish_cond) %>%PCA(quali.sup = c(1), quanti.sup = c(2), graph = FALSE)

dfish$pc1_allfish <- fish_pca$ind$coord[, 1]
dfish$pc2_allfish <- fish_pca$ind$coord[, 2]

pca.vars <- fish_pca$var$coord %>% data.frame
pca.vars$vars <- rownames(pca.vars)
pca.vars.m <- reshape2::melt(pca.vars, id.vars = "vars")
pca.var.sup <- fish_pca$quanti.sup$coord %>% data.frame

posArrowsX <- c(4, 4, 4, 4,4,4,4,4)
posArrowsY <- c(4, 4, 4, 4,4,4,4,4)
p_transf <- ggplot() +  
  geom_jitter(data = dfish, aes(x = pc1_allfish, y = pc2_allfish, colour = region_year), size = 2.5)+
  coord_equal()+
  theme_classic()+
  theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
        plot.background = element_rect(fill = "transparent", color = NA),
        legend.background = element_rect(fill=alpha('white', 0.4)))+
  labs(x= paste("PC1 (", round(fish_pca$eig["comp 1", "percentage of variance"], digits = 2), "%)"),
       y= paste("PC2 (", round(fish_pca$eig["comp 2", "percentage of variance"], digits = 2), "%)")) +
  geom_hline(yintercept = 0, lty = 2) +
  geom_vline(xintercept = 0, lty = 2) +
  geom_point(size = 2) +
  geom_segment(data = pca.vars, aes(x = 0, xend = Dim.1*4, y = 0, yend = Dim.2*4),
               arrow = arrow(length = unit(0.025, "npc"), type = "open"), 
               lwd = 0.5) + 
  geom_text(data = pca.vars, 
            aes(x = Dim.1*posArrowsX, y =  Dim.2*posArrowsY, label = vars), 
            check_overlap = FALSE, size = 5)+ 
  scale_color_manual(values = colours, name = "Region") +
  geom_segment(data = pca.var.sup, aes(x = 0, xend = Dim.1*5, y = 0, yend = Dim.2*5),
             arrow = arrow(length = unit(0.025, "npc"), type = "open"),
             lwd = 0.5, lty = "dotted") + 
  geom_text(data = pca.var.sup, 
            aes(x = Dim.1*4, y =  Dim.2*4,
                label = "Succès\nalimentaire"), size = 5, col = '#868B8D')
p_transf

ggsave(file="/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/PCA_fish.png", plot = p_transf, dpi = 500, width = 14, height = 10, bg = "transparent")

Les clusters se forment par station au minimum et au niveau de la région année, dans certains cas de la région-même. Donc la variation est vraiment due à l’environnement et aux poissons eux-mêmes plutôt qu’à la localisation elle-même. Et les clusters ne se distinguent pas entre l’Est et l’Ouest: tout est mélangé avec le Kitikmeot qui se distingue légèrement.

Contenus stomacaux

Cusa et al.

tot_preys2 <- prey_matrix_comb2 %>% select (c(1:11)) %>% rowSums()
prey_matrix_comb2[,1:11] <- prey_matrix_comb2[,1:11]/tot_preys2
  
FO_station <- prey_matrix_comb2 %>% aggregate(data = .,. ~ region_year+region+station, FUN = mean, na.rm = TRUE)

FO_station$shannon <- diversity(x = FO_station[,c(4:14)], index = "shannon", MARGIN = 1)
FO_station$simpson <- diversity(x = FO_station[,c(4:14)], index = "simpson", MARGIN = 1)

tot_carbon <- aggregate(carbon_mg ~ unique_fish_id + as.factor(station) + region_year, data = gutcontent2, FUN = sum)
colnames(tot_carbon) <- c("unique_fish_id", "station", "region_year", "sum_carbon")
A <- gutcontent2 %>% merge(tot_carbon, intersect(names(.), names(tot_carbon))) %>% 
  mutate(BDC_perc = carbon_mg/sum_carbon)

perc_carbon_fish_preys <- aggregate(BDC_perc ~ unique_fish_id + prey_id_comb + station + region_year + region, data= A, FUN = sum) %>% tidyr::spread(prey_id_comb, BDC_perc)
perc_carbon_fish_preys[is.na(perc_carbon_fish_preys)] <- 0

Xa_station <- perc_carbon_fish_preys %>% aggregate(data = .,. ~ region_year+station, FUN = mean, na.rm = TRUE)

Xap_station <- perc_carbon_fish_preys %>% aggregate(data = ., .~region_year+station, FUN = function(x){mean(x[x>0])})

is.nan.data.frame <- function(x){ do.call(cbind, lapply(x, is.nan))}
Xap_station[is.nan(Xap_station)] <- 0

env_reg_year <-  FO_station %>% select(1:2,15:22)%>%
  aggregate(. ~ region_year+region, data = ., FUN= mean, na.rm = TRUE)

preys_reg_year <- FO_station %>% select(-station)%>%
  aggregate(. ~ region_year+region, data = ., FUN= mean, na.rm = TRUE)
rownames(preys_reg_year)<- env_reg_year$region_year
rownames(env_reg_year)<- env_reg_year$region_year
carb_reg_year <- perc_carbon_fish_preys %>% aggregate(data = ., . ~ region_year+region, FUN = mean) %>% select(-unique_fish_id, -station)
rownames(carb_reg_year)<- env_reg_year$region_year

NMDS

MDS <- metaMDS(preys_reg_year[,3:13], try = 40)

## Run 0 stress 0.1592616 
## Run 1 stress 0.213972 
## Run 2 stress 0.1621729 
## Run 3 stress 0.2091287 
## Run 4 stress 0.1621731 
## Run 5 stress 0.1732711 
## Run 6 stress 0.3826912 
## Run 7 stress 0.2051073 
## Run 8 stress 0.2094691 
## Run 9 stress 0.2072162 
## Run 10 stress 0.2246623 
## Run 11 stress 0.1927286 
## Run 12 stress 0.1592616 
## ... New best solution
## ... Procrustes: rmse 3.360336e-05  max resid 8.203031e-05 
## ... Similar to previous best
## Run 13 stress 0.1732711 
## Run 14 stress 0.2114789 
## Run 15 stress 0.229833 
## Run 16 stress 0.1765413 
## Run 17 stress 0.1621736 
## Run 18 stress 0.1732717 
## Run 19 stress 0.2068415 
## Run 20 stress 0.1592618 
## ... Procrustes: rmse 0.0001301341  max resid 0.0004377828 
## ... Similar to previous best
## Run 21 stress 0.1621732 
## Run 22 stress 0.2213658 
## Run 23 stress 0.1592616 
## ... New best solution
## ... Procrustes: rmse 3.865589e-05  max resid 9.256398e-05 
## ... Similar to previous best
## Run 24 stress 0.2116018 
## Run 25 stress 0.1701745 
## Run 26 stress 0.1701744 
## Run 27 stress 0.1621727 
## Run 28 stress 0.1592616 
## ... New best solution
## ... Procrustes: rmse 1.437061e-05  max resid 3.994815e-05 
## ... Similar to previous best
## Run 29 stress 0.204023 
## Run 30 stress 0.1592629 
## ... Procrustes: rmse 0.0001749017  max resid 0.0004479572 
## ... Similar to previous best
## Run 31 stress 0.1621729 
## Run 32 stress 0.1621731 
## Run 33 stress 0.2225121 
## Run 34 stress 0.1592616 
## ... Procrustes: rmse 6.480072e-05  max resid 0.0002041358 
## ... Similar to previous best
## Run 35 stress 0.1732711 
## Run 36 stress 0.1621729 
## Run 37 stress 0.2177406 
## Run 38 stress 0.2267825 
## Run 39 stress 0.1732714 
## Run 40 stress 0.2059916 
## *** Solution reached

ord_env<- envfit(MDS ~ surf_temp_degC +surf_sal_kgm3 + open_water_day + NASC_zoo + indstandard + prof_mel, data=env_reg_year, perm=999)

par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
plot(MDS)
plot(ord_env)
orditorp(MDS,display="sites")
orditorp(MDS,display="species")
ordihull(MDS,
         group = preys_reg_year$region,
         col = c(1:8))

# legend("topright", inset = c(-0.8,0),
#        col = c(1:8), 
#        lty = 1,
#        legend = reg)

MDS2 <- metaMDS(carb_reg_year[,3:13], try = 40)

## Run 0 stress 0.1812863 
## Run 1 stress 0.1804913 
## ... New best solution
## ... Procrustes: rmse 0.1485507  max resid 0.5325009 
## Run 2 stress 0.1822983 
## Run 3 stress 0.1804914 
## ... Procrustes: rmse 5.225937e-05  max resid 0.0001488757 
## ... Similar to previous best
## Run 4 stress 0.1924862 
## Run 5 stress 0.1856083 
## Run 6 stress 0.2125887 
## Run 7 stress 0.1812863 
## Run 8 stress 0.1822983 
## Run 9 stress 0.1822986 
## Run 10 stress 0.2108407 
## Run 11 stress 0.1804921 
## ... Procrustes: rmse 0.001109986  max resid 0.004811274 
## ... Similar to previous best
## Run 12 stress 0.1804917 
## ... Procrustes: rmse 0.001018818  max resid 0.004463235 
## ... Similar to previous best
## Run 13 stress 0.1978196 
## Run 14 stress 0.1804942 
## ... Procrustes: rmse 0.001583977  max resid 0.006933139 
## ... Similar to previous best
## Run 15 stress 0.1804928 
## ... Procrustes: rmse 0.0004717767  max resid 0.002000221 
## ... Similar to previous best
## Run 16 stress 0.1822984 
## Run 17 stress 0.2175491 
## Run 18 stress 0.1822996 
## Run 19 stress 0.1822989 
## Run 20 stress 0.180491 
## ... New best solution
## ... Procrustes: rmse 0.000205268  max resid 0.0008695559 
## ... Similar to previous best
## Run 21 stress 0.1986603 
## Run 22 stress 0.1804928 
## ... Procrustes: rmse 0.0009966375  max resid 0.004324793 
## ... Similar to previous best
## Run 23 stress 0.1856082 
## Run 24 stress 0.1978179 
## Run 25 stress 0.1822985 
## Run 26 stress 0.1856083 
## Run 27 stress 0.1804909 
## ... New best solution
## ... Procrustes: rmse 9.642534e-05  max resid 0.0003826296 
## ... Similar to previous best
## Run 28 stress 0.1822983 
## Run 29 stress 0.1804941 
## ... Procrustes: rmse 0.001010252  max resid 0.004282965 
## ... Similar to previous best
## Run 30 stress 0.1822984 
## Run 31 stress 0.2340365 
## Run 32 stress 0.1804935 
## ... Procrustes: rmse 0.0008343145  max resid 0.00351976 
## ... Similar to previous best
## Run 33 stress 0.1804911 
## ... Procrustes: rmse 0.0001783757  max resid 0.0007399395 
## ... Similar to previous best
## Run 34 stress 0.1804938 
## ... Procrustes: rmse 0.0007932378  max resid 0.003311867 
## ... Similar to previous best
## Run 35 stress 0.1812863 
## Run 36 stress 0.1804911 
## ... Procrustes: rmse 0.000210149  max resid 0.0008782158 
## ... Similar to previous best
## Run 37 stress 0.1924919 
## Run 38 stress 0.1804909 
## ... Procrustes: rmse 5.51381e-05  max resid 0.00019667 
## ... Similar to previous best
## Run 39 stress 0.1804912 
## ... Procrustes: rmse 0.0002393642  max resid 0.0009951673 
## ... Similar to previous best
## Run 40 stress 0.2125813 
## *** Solution reached

ord_env2<- envfit(MDS2 ~ surf_temp_degC +surf_sal_kgm3 + open_water_day + NASC_zoo + indstandard + prof_mel, data=env_reg_year, perm=999)

par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
plot(MDS2)
plot(ord_env2)
orditorp(MDS2,display="sites", pch = 1)
orditorp(MDS2,display="species", pch = 17)
 ordihull(MDS2,
          group = carb_reg_year$region,
          col = c(1:8))

# legend("topright", inset = c(-0.8,0),
#        col = c(1:8), 
#        lty = 1,
#        legend = reg)

Clustering

bray_curt <- vegdist(preys_reg_year[,3:13], method = "bray")
plot(hclust(bray_curt), labels = unique(preys_reg_year$region_year))

adonis(bray_curt~ surf_temp_degC + surf_sal_kgm3 + open_water_day + NASC_zoo + indstandard + prof_mel, data = env_reg_year, strata = env_reg_year$region)

## 
## Call:
## adonis(formula = bray_curt ~ surf_temp_degC + surf_sal_kgm3 +      open_water_day + NASC_zoo + indstandard + prof_mel, data = env_reg_year,      strata = env_reg_year$region) 
## 
## Blocks:  strata 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)
## surf_temp_degC  1   0.15248 0.152477 1.75731 0.05751  0.405
## surf_sal_kgm3   1   0.17668 0.176680 2.03626 0.06664  0.699
## open_water_day  1   0.18528 0.185276 2.13532 0.06989  0.722
## NASC_zoo        1   0.11209 0.112090 1.29185 0.04228  0.291
## indstandard     1   0.12225 0.122249 1.40893 0.04611  0.882
## prof_mel        1   0.08022 0.080215 0.92449 0.03026  0.364
## Residuals      21   1.82211 0.086767         0.68730       
## Total          27   2.65109                  1.00000

adonis(bray_curt~ surf_temp_degC + surf_sal_kgm3 + open_water_day + NASC_zoo + indstandard + prof_mel, data = env_reg_year)

## 
## Call:
## adonis(formula = bray_curt ~ surf_temp_degC + surf_sal_kgm3 +      open_water_day + NASC_zoo + indstandard + prof_mel, data = env_reg_year) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)  
## surf_temp_degC  1   0.15248 0.152477 1.75731 0.05751  0.122  
## surf_sal_kgm3   1   0.17668 0.176680 2.03626 0.06664  0.078 .
## open_water_day  1   0.18528 0.185276 2.13532 0.06989  0.079 .
## NASC_zoo        1   0.11209 0.112090 1.29185 0.04228  0.260  
## indstandard     1   0.12225 0.122249 1.40893 0.04611  0.259  
## prof_mel        1   0.08022 0.080215 0.92449 0.03026  0.474  
## Residuals      21   1.82211 0.086767         0.68730         
## Total          27   2.65109                  1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

bray_curt_carb <- vegdist(carb_reg_year[,3:13], method = "bray")
plot(hclust(bray_curt_carb), labels = carb_reg_year$region_year)

adonis(bray_curt_carb ~ surf_temp_degC + surf_sal_kgm3 + open_water_day + NASC_zoo + indstandard + prof_mel, data = env_reg_year, strata = env_reg_year$region)

## 
## Call:
## adonis(formula = bray_curt_carb ~ surf_temp_degC + surf_sal_kgm3 +      open_water_day + NASC_zoo + indstandard + prof_mel, data = env_reg_year,      strata = env_reg_year$region) 
## 
## Blocks:  strata 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)  
## surf_temp_degC  1    0.2988 0.29883  3.0766 0.08421  0.370  
## surf_sal_kgm3   1    0.3036 0.30358  3.1255 0.08555  0.621  
## open_water_day  1    0.4044 0.40436  4.1630 0.11395  0.254  
## NASC_zoo        1    0.2691 0.26914  2.7709 0.07585  0.028 *
## indstandard     1    0.1193 0.11926  1.2278 0.03361  0.844  
## prof_mel        1    0.1135 0.11355  1.1690 0.03200  0.288  
## Residuals      21    2.0397 0.09713         0.57483         
## Total          27    3.5484                 1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis(bray_curt_carb ~ surf_temp_degC + surf_sal_kgm3 + open_water_day + NASC_zoo + indstandard + prof_mel, data = env_reg_year)

## 
## Call:
## adonis(formula = bray_curt_carb ~ surf_temp_degC + surf_sal_kgm3 +      open_water_day + NASC_zoo + indstandard + prof_mel, data = env_reg_year) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## surf_temp_degC  1    0.2988 0.29883  3.0766 0.08421  0.018 * 
## surf_sal_kgm3   1    0.3036 0.30358  3.1255 0.08555  0.010 **
## open_water_day  1    0.4044 0.40436  4.1630 0.11395  0.004 **
## NASC_zoo        1    0.2691 0.26914  2.7709 0.07585  0.017 * 
## indstandard     1    0.1193 0.11926  1.2278 0.03361  0.311   
## prof_mel        1    0.1135 0.11355  1.1690 0.03200  0.339   
## Residuals      21    2.0397 0.09713         0.57483          
## Total          27    3.5484                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

De ce que j’en conclus, le type de proies ingérées semble dépendre plus des conditions en tant que telles que la localisation géographique: les regions s’entremêlent dans l’espace. L’apport en carbone est modifié significativement par l’environnement, mais pas l’occurence des proies en tant que tel.

Section 2. Modèles de succès alimentaire

Explo préliminaire

#taille des poissons
plot(est_standard_length ~ region, data = dfish, las = 2)

tapply(dfish$est_standard_length, INDEX = list(dfish$region, dfish$year), FUN = mean)

##                                   2009     2010     2011     2014     2015
## Amundsen Gulf Mouth                 NA       NA       NA       NA 32.85720
## Coronation Maud                     NA       NA 15.59700       NA       NA
## Lancaster Sound                     NA 13.23085 13.62750       NA 16.59357
## Larsen Sound - Victoria Strait      NA 11.30658       NA       NA       NA
## Mackenzie Shelf                12.6857 27.09488       NA 35.40000 31.54290
## NEG                                 NA       NA       NA       NA       NA
## North Water                         NA       NA       NA 17.73251       NA
## Peel Sound                          NA 10.81823 17.42556 12.13674 11.36489
## West Baffin Bay                     NA       NA       NA       NA 14.48783
##                                    2016     2017     2018
## Amundsen Gulf Mouth                  NA       NA       NA
## Coronation Maud                22.55682 16.26570 14.98442
## Lancaster Sound                17.96064 15.33450       NA
## Larsen Sound - Victoria Strait 14.74777       NA 11.03090
## Mackenzie Shelf                      NA       NA       NA
## NEG                                  NA 19.12333       NA
## North Water                    20.41569       NA 21.30636
## Peel Sound                     12.85020       NA       NA
## West Baffin Bay                11.19310       NA       NA

tapply(dfish$est_standard_length, INDEX = list(dfish$region, dfish$year), FUN = median)

##                                  2009   2010  2011   2014   2015   2016    2017
## Amundsen Gulf Mouth                NA     NA    NA     NA 35.000     NA      NA
## Coronation Maud                    NA     NA 16.07     NA     NA 18.906 16.1720
## Lancaster Sound                    NA 12.713 13.14     NA 15.857 18.438 14.6875
## Larsen Sound - Victoria Strait     NA 11.431    NA     NA     NA 14.063      NA
## Mackenzie Shelf                12.307 27.259    NA 38.000 33.500     NA      NA
## NEG                                NA     NA    NA     NA     NA     NA 18.6500
## North Water                        NA     NA    NA 18.924     NA 21.997      NA
## Peel Sound                         NA 10.790 12.14 12.118 11.286 13.047      NA
## West Baffin Bay                    NA     NA    NA     NA 13.857 11.136      NA
##                                   2018
## Amundsen Gulf Mouth                 NA
## Coronation Maud                15.1560
## Lancaster Sound                     NA
## Larsen Sound - Victoria Strait 10.7025
## Mackenzie Shelf                     NA
## NEG                                 NA
## North Water                    21.8750
## Peel Sound                          NA
## West Baffin Bay                     NA

cor(dfish$est_standard_length, dfish$open_water_day)

## [1] 0.4979243

plot(dfish$est_standard_length~dfish$open_water_day)

lattice::xyplot(est_standard_length ~ open_water_day|region, data = dfish)

#taille et carbone
plot(dfish$carbon_mg ~ dfish$fish_WW)

plot(dfish$carbon_mg ~ dfish$est_standard_length) #la quantité de carbone ingérée demeure faible même si les poissons sont gros, mais certains gros individus vont arriver à toucher le jackpot

#taille et succès alimentaire
plot((feeding_success)~est_standard_length, data = dfish)

plot(log10(feeding_success)~est_standard_length, data = dfish)

summary(lm(log10(feeding_success) ~ est_standard_length, data = dfish)) # il va falloir corriger pour la taille dans les modèles: elle semble influencer le succès alimentaire

## 
## Call:
## lm(formula = log10(feeding_success) ~ est_standard_length, data = dfish)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.42090 -0.26394  0.05922  0.30865  1.15647 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -3.111077   0.058619 -53.073  < 2e-16 ***
## est_standard_length  0.009872   0.003030   3.258  0.00124 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4481 on 337 degrees of freedom
## Multiple R-squared:  0.03053,    Adjusted R-squared:  0.02766 
## F-statistic: 10.61 on 1 and 337 DF,  p-value: 0.001237

plot(log10(feeding_success) ~ est_standard_length, data= dfish[dfish$est_standard_length>=25,]) 

plot(log10(feeding_success) ~ est_standard_length, data= dfish[dfish$est_standard_length<25,])

plot(fish_cond ~ est_standard_length,data = dfish) # les poissons plus gros n'ont pas nécessairement une meilleure condition physique, et je ne discerne pas vraiment un pattern indiquant une transition vers une autre forme

plot(log10(feeding_success) ~ fish_cond, data  = dfish) #augmentation de FS

plot(log10(feeding_success) ~ fish_cond, data= dfish[dfish$est_standard_length>=25,])

plot(log10(feeding_success) ~ fish_cond, data= dfish[dfish$est_standard_length<25,]) 

plot(fish_cond~feeding_success, data = dfish)

plot(fish_cond~log10(feeding_success), data = dfish) # légère augmentation quand FS augmente

dfish$size_catego <- with(dfish, ifelse(est_standard_length<10, "<10", 
                                        ifelse(est_standard_length >=10 & est_standard_length<15, "10-15",
                                               ifelse(est_standard_length>=15 & est_standard_length<20, "15-20",
                                                      ifelse(est_standard_length>=20 & est_standard_length<25, "20-25",
                                                             ">25")))))
dfish$logFS <- log10(dfish$feeding_success)
xyplot(logFS ~ fish_cond|size_catego, data = dfish)

#FS et environnement
plot(logFS ~ NASC_zoo, data = dfish) #très légère augmentation mais rien de signifcatif

plot(logFS ~ surf_temp_degC, data = dfish) #un peu en parabole

plot(logFS ~ surf_sal_kgm3, data = dfish) #légère augmentation

plot(logFS ~ open_water_day, data = dfish) #augmentation

plot(logFS ~ indstandard, data = dfish)

dfish %>% select(logFS, open_water_day, surf_sal_kgm3, surf_temp_degC, NASC_zoo, fish_cond, indstandard) %>% plot

dfish %>% select(logFS, open_water_day, surf_sal_kgm3, surf_temp_degC, NASC_zoo, fish_cond, indstandard) %>% cor

##                       logFS open_water_day surf_sal_kgm3 surf_temp_degC
## logFS           1.000000000     0.28799513    0.19917441     0.04315438
## open_water_day  0.287995134     1.00000000    0.46112232     0.35527157
## surf_sal_kgm3   0.199174408     0.46112232    1.00000000    -0.16953030
## surf_temp_degC  0.043154376     0.35527157   -0.16953030     1.00000000
## NASC_zoo       -0.003297709     0.45518277    0.24197166     0.56389798
## fish_cond       0.170762651     0.15052952   -0.08272398     0.16134722
## indstandard     0.042107243     0.07052434    0.13824947    -0.21632976
##                    NASC_zoo   fish_cond indstandard
## logFS          -0.003297709  0.17076265  0.04210724
## open_water_day  0.455182768  0.15052952  0.07052434
## surf_sal_kgm3   0.241971663 -0.08272398  0.13824947
## surf_temp_degC  0.563897983  0.16134722 -0.21632976
## NASC_zoo        1.000000000 -0.10299265 -0.15853361
## fish_cond      -0.102992654  1.00000000 -0.11306133
## indstandard    -0.158533609 -0.11306133  1.00000000

#corrélation entre  débâcle et NASC_zoo, NASC_zoo et température
boxplot(logFS ~ region,data = dfish)

boxplot(feeding_success ~ region_year,data = dfish, las = 2)

boxplot(logFS ~ region_year,data = dfish, las = 2)

boxplot(logFS ~ province, data = dfish)

plot(feeding_success ~ region, data = dfish, las = 2) #pas mal similaire, mais toutes des valeurs bound à 0

plot(log10(feeding_success) ~ region, data = dfish, las = 2)

hist(dfish$feeding_success) #un beau neg binomial

hist(log10(dfish$feeding_success)) #normalish

lattice::bwplot(feeding_success ~ as.factor(year)|region, data = dfish)

lattice::bwplot(log10(feeding_success) ~ as.factor(year)|region, data = dfish) # variation LS, LV, MS, NW, PS

lattice::xyplot(log10(feeding_success) ~ (open_water_day)|region, data = dfish) #la variation individuelle est énorme toute région confondue

plot(feeding_success ~ as.factor(province), data = dfish)

plot(log10(feeding_success) ~ as.factor(province), data = dfish)

FS ~ environnement

#mixtes linéraires avec transfo
#modèle zoo et compétition
zoo_fish <- lme(logFS~indstandard + NASC_zoo, random = ~1|region, data = dfish)
summary(zoo_fish) #la compétition et l'abondance de macrozooplancton n'ont pas d'effet sur le succès alimentaire

## Linear mixed-effects model fit by REML
##  Data: dfish 
##        AIC     BIC    logLik
##   432.0875 451.173 -211.0437
## 
## Random effects:
##  Formula: ~1 | region
##         (Intercept)  Residual
## StdDev:   0.1646009 0.4317488
## 
## Fixed effects: logFS ~ indstandard + NASC_zoo 
##                  Value  Std.Error  DF   t-value p-value
## (Intercept) -2.8918384 0.06779990 328 -42.65254  0.0000
## indstandard  0.0368243 0.06684152 328   0.55092  0.5821
## NASC_zoo    -0.0043669 0.00392956 328  -1.11130  0.2673
##  Correlation: 
##             (Intr) indstn
## indstandard -0.261       
## NASC_zoo    -0.374  0.102
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.2312115 -0.5861180  0.1224188  0.6502431  3.0791358 
## 
## Number of Observations: 339
## Number of Groups: 9

#modèle débâcle et température
bu_temp <- lme(logFS~surf_temp_degC + surf_sal_kgm3 + open_water_day, random = ~1|region, data = dfish)
summary(bu_temp) # la salinité et la date de débâcle influencent le succès alimentaire

## Linear mixed-effects model fit by REML
##  Data: dfish 
##        AIC      BIC    logLik
##   420.2667 443.1515 -204.1333
## 
## Random effects:
##  Formula: ~1 | region
##         (Intercept)  Residual
## StdDev:   0.2367474 0.4124639
## 
## Fixed effects: logFS ~ surf_temp_degC + surf_sal_kgm3 + open_water_day 
##                    Value Std.Error  DF    t-value p-value
## (Intercept)    -3.944336 0.3376341 327 -11.682280  0.0000
## surf_temp_degC -0.022784 0.0136453 327  -1.669716  0.0959
## surf_sal_kgm3   0.025262 0.0109628 327   2.304337  0.0218
## open_water_day  0.007498 0.0015022 327   4.991484  0.0000
##  Correlation: 
##                (Intr) srf__C srf__3
## surf_temp_degC -0.237              
## surf_sal_kgm3  -0.943  0.210       
## open_water_day -0.168 -0.344 -0.024
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.1991235 -0.5861290  0.0786991  0.6872109  3.1684312 
## 
## Number of Observations: 339
## Number of Groups: 9

#modèle de taille et condition physique
size_cond <- lme(logFS~est_standard_length + fish_cond, random = ~1|region, data = dfish)
summary(size_cond) #la taille et la condition physique influencent le succès alimentaire

## Linear mixed-effects model fit by REML
##  Data: dfish 
##        AIC      BIC    logLik
##   414.5845 433.6701 -202.2923
## 
## Random effects:
##  Formula: ~1 | region
##         (Intercept)  Residual
## StdDev:   0.1495965 0.4212912
## 
## Fixed effects: logFS ~ est_standard_length + fish_cond 
##                          Value  Std.Error  DF   t-value p-value
## (Intercept)         -3.1493222 0.08503158 328 -37.03709  0.0000
## est_standard_length  0.0129550 0.00348403 328   3.71840  0.0002
## fish_cond            0.1983006 0.07598673 328   2.60967  0.0095
##  Correlation: 
##                     (Intr) est_s_
## est_standard_length -0.75        
## fish_cond           -0.05   0.06 
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.9914524 -0.5812885  0.1124799  0.6265541  3.1769199 
## 
## Number of Observations: 339
## Number of Groups: 9

FS_cond <- lme(fish_cond~feeding_success, random = ~1|region, data = dfish)
summary(FS_cond) #le succès alimentaire n'a pas d'effet sur la condition physique

## Linear mixed-effects model fit by REML
##  Data: dfish 
##        AIC      BIC    logLik
##   165.8988 181.1792 -78.94942
## 
## Random effects:
##  Formula: ~1 | region
##         (Intercept)  Residual
## StdDev:   0.1559522 0.2980211
## 
## Fixed effects: fish_cond ~ feeding_success 
##                     Value Std.Error  DF    t-value p-value
## (Intercept)     -0.019643  0.057479 329 -0.3417334  0.7328
## feeding_success 14.865890  8.345654 329  1.7812731  0.0758
##  Correlation: 
##                 (Intr)
## feeding_success -0.274
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.582522196 -0.629106889 -0.008520807  0.588758602  4.156256376 
## 
## Number of Observations: 339
## Number of Groups: 9

#correction pour taille
bu_temp2 <- lme(logFS~surf_temp_degC + surf_sal_kgm3 + open_water_day + fish_cond, random = list(~1|region, ~1|size_catego), data = dfish)
summary(bu_temp2) # salinité, condition physique et date de débâcle ont un effet sur le succès alimentaire, même lorsqu'on corrige pour la taille

## Linear mixed-effects model fit by REML
##  Data: dfish 
##        AIC      BIC   logLik
##   422.2459 452.7351 -203.123
## 
## Random effects:
##  Formula: ~1 | region
##         (Intercept)
## StdDev:    0.226057
## 
##  Formula: ~1 | size_catego %in% region
##         (Intercept)  Residual
## StdDev:  0.05983658 0.4077269
## 
## Fixed effects: logFS ~ surf_temp_degC + surf_sal_kgm3 + open_water_day + fish_cond 
##                    Value Std.Error  DF    t-value p-value
## (Intercept)    -3.930930 0.3363123 299 -11.688332  0.0000
## surf_temp_degC -0.023752 0.0135917 299  -1.747532  0.0816
## surf_sal_kgm3   0.025785 0.0109281 299   2.359504  0.0189
## open_water_day  0.007114 0.0014886 299   4.779315  0.0000
## fish_cond       0.170881 0.0751147 299   2.274936  0.0236
##  Correlation: 
##                (Intr) srf__C srf__3 opn_w_
## surf_temp_degC -0.249                     
## surf_sal_kgm3  -0.945  0.218              
## open_water_day -0.165 -0.331 -0.027       
## fish_cond      -0.031 -0.025  0.049 -0.055
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.24599508 -0.59573689  0.07092356  0.69214517  3.27656891 
## 
## Number of Observations: 339
## Number of Groups: 
##                  region size_catego %in% region 
##                       9                      36

zoo_fish2 <- lme(logFS~indstandard + NASC_zoo, random = list(~1|region, ~1|size_catego), data = dfish)
summary(zoo_fish2) #pas d'effet de la compétition ni de l'abondance de macrozooplancton

## Linear mixed-effects model fit by REML
##  Data: dfish 
##        AIC      BIC    logLik
##   431.9769 454.8796 -209.9885
## 
## Random effects:
##  Formula: ~1 | region
##         (Intercept)
## StdDev:   0.1470768
## 
##  Formula: ~1 | size_catego %in% region
##         (Intercept)  Residual
## StdDev:  0.09833968 0.4253594
## 
## Fixed effects: logFS ~ indstandard + NASC_zoo 
##                  Value  Std.Error  DF   t-value p-value
## (Intercept) -2.8707589 0.06643741 301 -43.20998  0.0000
## indstandard  0.0329818 0.06555735 301   0.50310  0.6153
## NASC_zoo    -0.0052195 0.00388055 301  -1.34504  0.1796
##  Correlation: 
##             (Intr) indstn
## indstandard -0.264       
## NASC_zoo    -0.390  0.105
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.3095903 -0.5411252  0.1185393  0.6376297  3.1496849 
## 
## Number of Observations: 339
## Number of Groups: 
##                  region size_catego %in% region 
##                       9                      36

Section 3. Proies et succès alimentaire

FS ~ taille des proies

plot(prey_range_um/1000 ~ est_standard_length, data = dfish) # les grosses larves mangent à la fois de grosses et petites choses

plot(prey_min_um/1000 ~ est_standard_length, data = dfish)

plot(prey_max_um/1000 ~ est_standard_length, data = dfish)

plot(prey_median_um/1000 ~ est_standard_length, data = dfish) #vraiment rare que les larves, peu importe la taille, mangent de grosses affaires

plot(log10(feeding_success) ~ prey_median_um, data = dfish)

median(dfish[which(dfish$feeding_success > 0.001),"prey_median_um"])

## [1] 269.5

plot(log10(feeding_success) ~ prey_max_um, data = dfish)

plot(log10(feeding_success) ~ prey_min_um, data = dfish)

plot(log10(feeding_success) ~ prey_range_um, data = dfish)

plot(log10(feeding_success) ~ prey_median_um, data = dfish)

plot((feeding_success) ~ prey_max_um, data = dfish)

plot((feeding_success) ~ prey_min_um, data = dfish)

plot((feeding_success) ~ prey_range_um, data = dfish)

plot((feeding_success) ~ prey_median_um, data = dfish)

#pas vraiment de relation entre taille des proies et FS

dfish$total_preys <- rowSums(dfish[,c(23:86)])
plot(total_preys ~ prey_max_um, data = dfish)

plot(total_preys ~ prey_min_um, data = dfish)

plot(total_preys~prey_range_um, data = dfish)

plot(total_preys~prey_median_um, data = dfish)

plot((feeding_success) ~ total_preys, data = dfish)

plot(log10(feeding_success) ~ total_preys, data = dfish)

FS_preySize <- lme(logFS ~ prey_median_um + prey_range_um + total_preys, random = list(~1|region, ~1|size_catego), data = dfish)
summary(FS_preySize)

## Linear mixed-effects model fit by REML
##  Data: dfish 
##        AIC      BIC    logLik
##   444.1923 470.8912 -215.0962
## 
## Random effects:
##  Formula: ~1 | region
##         (Intercept)
## StdDev:   0.1348462
## 
##  Formula: ~1 | size_catego %in% region
##         (Intercept)  Residual
## StdDev:  0.05866081 0.4124322
## 
## Fixed effects: logFS ~ prey_median_um + prey_range_um + total_preys 
##                     Value  Std.Error  DF   t-value p-value
## (Intercept)    -2.9160423 0.06843441 300 -42.61076  0.0000
## prey_median_um -0.0002129 0.00010534 300  -2.02075  0.0442
## prey_range_um  -0.0000138 0.00002629 300  -0.52648  0.5989
## total_preys     0.0011068 0.00023005 300   4.81124  0.0000
##  Correlation: 
##                (Intr) pry_m_ pry_r_
## prey_median_um -0.419              
## prey_range_um  -0.169 -0.388       
## total_preys    -0.347  0.061  0.024
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.16006905 -0.56156550  0.08504287  0.64312527  2.96547955 
## 
## Number of Observations: 339
## Number of Groups: 
##                  region size_catego %in% region 
##                       9                      36

FS ~ prés absence de proies

prey_matrix_comb3 <- prey_matrix_comb2
prey_matrix_comb3[, c(1:11)] <- as.data.frame(sapply(prey_matrix_comb3[,c(1:11)], as.logical))
prey_matrix_comb3[, c(1:11)] <- as.data.frame(sapply(prey_matrix_comb3[,c(1:11)], as.numeric))

plot(prey_matrix_comb3[, c(1:11, 23)])

summary(lm(log10(feeding_success) ~ appendicularia + `calanus glacialis c` +`calanus sp n`+`cyclopoidae sp c`+`egg`+`other` + `other calanoid sp c`+`other calanus sp c`+`other copepodite and n`+`pseudocalanus sp c`+`pseudocalanus sp n`, data = prey_matrix_comb3))

## 
## Call:
## lm(formula = log10(feeding_success) ~ appendicularia + `calanus glacialis c` + 
##     `calanus sp n` + `cyclopoidae sp c` + egg + other + `other calanoid sp c` + 
##     `other calanus sp c` + `other copepodite and n` + `pseudocalanus sp c` + 
##     `pseudocalanus sp n`, data = prey_matrix_comb3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.97309 -0.22962  0.00821  0.22896  1.37086 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              -3.432620   0.070979 -48.361  < 2e-16 ***
## appendicularia            0.142919   0.075528   1.892  0.05934 .  
## `calanus glacialis c`     0.149581   0.052713   2.838  0.00483 ** 
## `calanus sp n`            0.349190   0.058711   5.948 6.99e-09 ***
## `cyclopoidae sp c`       -0.109163   0.044016  -2.480  0.01364 *  
## egg                      -0.001433   0.050067  -0.029  0.97718    
## other                    -0.016998   0.044173  -0.385  0.70063    
## `other calanoid sp c`     0.033220   0.048880   0.680  0.49722    
## `other calanus sp c`      0.292917   0.051564   5.681 2.97e-08 ***
## `other copepodite and n` -0.056865   0.043665  -1.302  0.19373    
## `pseudocalanus sp c`      0.213253   0.046922   4.545 7.75e-06 ***
## `pseudocalanus sp n`      0.053544   0.052531   1.019  0.30882    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3645 on 327 degrees of freedom
## Multiple R-squared:  0.3777, Adjusted R-squared:  0.3568 
## F-statistic: 18.04 on 11 and 327 DF,  p-value: < 2.2e-16

colnames(prey_matrix_comb3) <- make.names(colnames(prey_matrix_comb3), unique = TRUE)
prey_matrix_comb3$logFS <- log10(prey_matrix_comb3$feeding_success)
summary(lme(as.formula(paste(colnames(prey_matrix_comb3)[24], "~",
        paste(colnames(prey_matrix_comb3)[c(1:11)], collapse = "+"),
        sep = "")), random = ~1|region, data = prey_matrix_comb3))

## Linear mixed-effects model fit by REML
##  Data: prey_matrix_comb3 
##       AIC      BIC    logLik
##   326.237 379.2964 -149.1185
## 
## Random effects:
##  Formula: ~1 | region
##         (Intercept)  Residual
## StdDev:   0.1262702 0.3468517
## 
## Fixed effects: as.formula(paste(colnames(prey_matrix_comb3)[24], "~", paste(colnames(prey_matrix_comb3)[c(1:11)],      collapse = "+"), sep = "")) 
##                            Value  Std.Error  DF   t-value p-value
## (Intercept)            -3.423451 0.08259986 319 -41.44621  0.0000
## appendicularia          0.154592 0.07352840 319   2.10248  0.0363
## calanus.glacialis.c     0.132308 0.05273530 319   2.50890  0.0126
## calanus.sp.n            0.325550 0.05611091 319   5.80190  0.0000
## cyclopoidae.sp.c       -0.091938 0.04362995 319  -2.10722  0.0359
## egg                     0.010304 0.04855142 319   0.21222  0.8321
## other                  -0.041554 0.04385643 319  -0.94750  0.3441
## other.calanoid.sp.c     0.037001 0.04716550 319   0.78448  0.4333
## other.calanus.sp.c      0.299861 0.05103932 319   5.87510  0.0000
## other.copepodite.and.n -0.005662 0.04304434 319  -0.13153  0.8954
## pseudocalanus.sp.c      0.224864 0.04643036 319   4.84303  0.0000
## pseudocalanus.sp.n      0.017630 0.05161680 319   0.34155  0.7329
##  Correlation: 
##                        (Intr) appndc clns.g. clns.s. cycl.. egg    other 
## appendicularia          0.058                                            
## calanus.glacialis.c    -0.106 -0.068                                     
## calanus.sp.n           -0.344 -0.023 -0.015                              
## cyclopoidae.sp.c       -0.009 -0.054 -0.021  -0.080                      
## egg                    -0.388 -0.035  0.037  -0.045  -0.058              
## other                  -0.175 -0.099 -0.037   0.088  -0.233 -0.136       
## other.calanoid.sp.c    -0.021 -0.030 -0.100  -0.057  -0.028  0.101 -0.158
## other.calanus.sp.c     -0.144 -0.038 -0.363   0.045  -0.023  0.078  0.070
## other.copepodite.and.n -0.152 -0.105  0.109  -0.113  -0.096 -0.042 -0.036
## pseudocalanus.sp.c     -0.089  0.067 -0.237  -0.109  -0.146  0.068 -0.039
## pseudocalanus.sp.n     -0.245 -0.078  0.095  -0.261  -0.137 -0.072  0.122
##                        othr.clnd.. othr.clns.. othr.cp.. psdclns.sp.c
## appendicularia                                                       
## calanus.glacialis.c                                                  
## calanus.sp.n                                                         
## cyclopoidae.sp.c                                                     
## egg                                                                  
## other                                                                
## other.calanoid.sp.c                                                  
## other.calanus.sp.c     -0.152                                        
## other.copepodite.and.n -0.005      -0.035                            
## pseudocalanus.sp.c     -0.238      -0.071      -0.019                
## pseudocalanus.sp.n     -0.017       0.046      -0.137    -0.081      
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.36210098 -0.58662279  0.01303377  0.62284450  3.57345605 
## 
## Number of Observations: 339
## Number of Groups: 9

prey_matrix_comb3$size_catego <- with(prey_matrix_comb3, ifelse(est_standard_length<10, "<10", 
                                        ifelse(est_standard_length >=10 & est_standard_length<15, "10-15",
                                               ifelse(est_standard_length>=15 & est_standard_length<20, "15-20",
                                                      ifelse(est_standard_length>=20 & est_standard_length<25, "20-25",
                                                             ">25")))))
FS_preyType <- lme(as.formula(paste(colnames(prey_matrix_comb3)[24], "~",
        paste(colnames(prey_matrix_comb3)[c(1:11)], collapse = "+"),
        sep = "")), random = c(~1|region, ~1|size_catego), data = prey_matrix_comb3)
summary(FS_preyType)

## Linear mixed-effects model fit by REML
##  Data: prey_matrix_comb3 
##        AIC      BIC    logLik
##   327.7307 384.5801 -148.8654
## 
## Random effects:
##  Formula: ~1 | region
##         (Intercept)
## StdDev:   0.1161706
## 
##  Formula: ~1 | size_catego %in% region
##         (Intercept)  Residual
## StdDev:  0.05383227 0.3447138
## 
## Fixed effects: as.formula(paste(colnames(prey_matrix_comb3)[24], "~", paste(colnames(prey_matrix_comb3)[c(1:11)],      collapse = "+"), sep = "")) 
##                            Value  Std.Error  DF   t-value p-value
## (Intercept)            -3.427844 0.08193596 292 -41.83565  0.0000
## appendicularia          0.152070 0.07343947 292   2.07069  0.0393
## calanus.glacialis.c     0.137943 0.05320298 292   2.59276  0.0100
## calanus.sp.n            0.327780 0.05595309 292   5.85813  0.0000
## cyclopoidae.sp.c       -0.089545 0.04364232 292  -2.05179  0.0411
## egg                     0.014573 0.04851072 292   0.30041  0.7641
## other                  -0.038711 0.04388143 292  -0.88216  0.3784
## other.calanoid.sp.c     0.040297 0.04710597 292   0.85546  0.3930
## other.calanus.sp.c      0.302086 0.05102480 292   5.92038  0.0000
## other.copepodite.and.n -0.006180 0.04289898 292  -0.14407  0.8855
## pseudocalanus.sp.c      0.223692 0.04645031 292   4.81572  0.0000
## pseudocalanus.sp.n      0.014466 0.05169749 292   0.27982  0.7798
##  Correlation: 
##                        (Intr) appndc clns.g. clns.s. cycl.. egg    other 
## appendicularia          0.055                                            
## calanus.glacialis.c    -0.107 -0.063                                     
## calanus.sp.n           -0.347 -0.022 -0.014                              
## cyclopoidae.sp.c       -0.013 -0.051 -0.017  -0.080                      
## egg                    -0.390 -0.035  0.037  -0.046  -0.053              
## other                  -0.179 -0.098 -0.036   0.086  -0.231 -0.137       
## other.calanoid.sp.c    -0.027 -0.028 -0.093  -0.056  -0.021  0.103 -0.158
## other.calanus.sp.c     -0.149 -0.033 -0.349   0.048  -0.020  0.076  0.075
## other.copepodite.and.n -0.151 -0.104  0.109  -0.112  -0.100 -0.044 -0.031
## pseudocalanus.sp.c     -0.094  0.065 -0.245  -0.105  -0.144  0.066 -0.035
## pseudocalanus.sp.n     -0.246 -0.077  0.079  -0.262  -0.136 -0.072  0.120
##                        othr.clnd.. othr.clns.. othr.cp.. psdclns.sp.c
## appendicularia                                                       
## calanus.glacialis.c                                                  
## calanus.sp.n                                                         
## cyclopoidae.sp.c                                                     
## egg                                                                  
## other                                                                
## other.calanoid.sp.c                                                  
## other.calanus.sp.c     -0.147                                        
## other.copepodite.and.n -0.005      -0.032                            
## pseudocalanus.sp.c     -0.236      -0.070      -0.022                
## pseudocalanus.sp.n     -0.019       0.040      -0.137    -0.076      
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -2.489743779 -0.568567330  0.008679755  0.623507919  3.534221547 
## 
## Number of Observations: 339
## Number of Groups: 
##                  region size_catego %in% region 
##                       9                      36